# Implement the rescaled change time method on single-day data
# For AAPL 2019-01-02, WITHSIZE, SUPPORT 4, DELTA 0.5

library(data.table)
library(tidyr)
library(xts)
library(dplyr)
library(pbmcapply)
load("~/Desktop/orderbook/Research/R/gof/AAPL_2019-01-02.RData")

support_max = 20
delta = 0.25
config = "nosize"
#config = ""
if(config == "nosize"){
  event_size = ""
}else{
  event_size = "eventsize_"
}

model_config = "_state_hawkes_time_customizedtwo_LASSO_0005"
# Load Hawkes estimation data
hawkes = readRDS(paste0("~/Desktop/orderbook/Research/R/gof/images/hawkes_markov_pref_price_", event_size,
                        data_file_name, paste0("_support",toString(support_max),"delta", toString(delta),model_config,".rds") ))


delta = hawkes$delta
support_max = hawkes$support_max
support_seq = seq(0,support_max,delta)

book =  book[modifying_queue %in% -3:3]

book = book[, c("time","queue_and_event")]


# Selected event for goodness-of-fit testing is "+2(+)"
event_names = colnames(event)[-1]


pdf(paste0("~/Desktop/Research/R/gof/images/gof_plots_all_events_", data_file_name, ".pdf"), width = 10, height = 5    )

for (i in 1:1){
  selected_event = event_names[i]
  selected_book = book[queue_and_event == selected_event]
  
  empirical_1 = readRDS(paste0("~/Desktop/orderbook/Research/R/gof/images/empirical_",event_size,data_file_name,"_support",toString(support_max), "delta", "0.5",model_config,"_",selected_event,".rds" ))
  empirical_2 = readRDS(paste0("~/Desktop/orderbook/Research/R/gof/images/empirical_",event_size,data_file_name,"_support",toString(support_max), "delta", "0.125",model_config,"_",selected_event,".rds" ))
  empirical_3 = readRDS(paste0("~/Desktop/orderbook/Research/R/gof/images/empirical_",event_size,data_file_name,"_support",toString(support_max), "delta", "0.5","_hawkes_LASSO_0005","_",selected_event,".rds" ))
  
  print(mean(empirical_1))
  print(var(empirical_1))
  
  print(mean(empirical_2))
  print(var(empirical_2))
  
  print(mean(empirical_3))
  print(var(empirical_3))
  
  homo_intensity = dim(selected_book)[1] / (max(book$time) - min(book$time))
  empirical_homo = diff(selected_book$time) * homo_intensity
  print(mean(empirical_homo))
  print(var(empirical_homo))
  
  
  p <- ppoints(100)    # 100 equally spaced points on (0,1), excluding endpoints
  q_homo <- quantile(empirical_homo,p=p) # percentiles of the sample distribution
  q1 = quantile(empirical_1, p=p)
  q2 = quantile(empirical_2, p=p)
  q3 = quantile(empirical_3, p=p)
  q_std_exp = quantile( rexp(100000),p=p    )
  
  par(oma=c(0, 0, 0, 10))
  plot(log10(qexp(p)) , log10(q_homo), main=paste0("event ",selected_event),
       xlab="Theoretical Exponential Quantiles",ylab=paste0("LOB Data Quantiles for Event ", selected_event),cex = 0.2)
  points(log10(qexp(p)),log10(q1),col="red" ,cex = 0.2)
  points(log10(qexp(p)),log10(q2),col="green" ,cex = 0.2)
  points(log10(qexp(p)),log10(q3),col="purple" ,cex = 0.2)
  points(log10(qexp(p)),log10(q_std_exp),col="blue", cex= 0.2 )
  
  legend(par('usr')[2], 
         bty = "n",
         legend=c("Homogeneous Poisson","Delta 0.5s","Delta 0.125s","Delta 0.5s hawkes only","Unit exponential dist"), 
         pch=c(1,1,1,1,1), 
         col =c("black","red","green","purple","blue"),
         xpd =NA,cex= 0.7)
}
dev.off()
